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Abstract 



We re-examine the procedure of adiabatic elimination of fast relaxing vari- 
ables near a bifurcation point when some of the parameters of the system are 
stochastically modulated. Approximate stationary solutions of the Fokker- 
Planck equation are obtained near threshold for the pitchfork and transcritical 
bifurcations. Stochastic resonance between fast variables and random mod- 
ulation may shift the effective bifurcation point by an amount proportional 
to the intensity of the fluctuations. We also find that fluctuations of the 
fast variables above threshold are not always Gaussian and centered around 
the (deterministic) center manifold as was previously believed. Numerical 
solutions obtained for a few illustrative examples support these conclusions. 
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Typeset using REVT^ 



I. INTRODUCTION 



A system is said to undergo a bifurcation when its long time behavior changes quahta- 
tively as some control parameter is continuously varied. Examples include the saddle-node, 
transcritical and pitchfork bifurcations, which involve a transition between two fixed point 
solutions, and the Hopf bifurcation that involves a transition between a fixed point solution 
and a limit cycle. Near the bifurcation point only a small number of so-called slow variables 
are required to determine the evolution of the system over a long time scale. The remaining 
degrees of freedom (the so-called fast variables) adjust very rapidly to the instantaneous 
values of the slow variables, and can be adiabatically eliminated. The qualitative features of 
the evolution of the system near the bifurcation point are thus obtained by constraining the 
original governing equations to a surface in phase space known as the center manifold. The 
resulting equations valid on the manifold are the normal form equations The purpose 
of this article is to re-examine the analogous reduction procedure when one or more of the 
system's parameters include a random component [0-0 . 

We focus mainly on the case in which the externally set control parameter includes a 
small random component which we model as a stochastic process in time. In this case, the 
bifurcation point remains sharp, although its location may depend on the intensity of the 
fluctuations. Although there is arguably little conceptual difference between deterministic 
variables that relax quickly in the vicinity of the bifurcation point, and a stochastic process 
of short correlation time (say of the same order or smaller than inverse relaxation rates of 
the fast variables), we show below that stochastic resonance between the two can affect the 
evolution on the slow time scale. 

The essential aspects of the adiabatic reduction procedure in the stochastic case can 
be illustrated in the simple case of a second order system. Let A be the amplitude of a 
bifurcating mode, and B the amplitude of a second mode that is itself linearly stable near 
onset. A reduced control parameter a is defined such that the trivial state A = B = is 
stable if a < 0, and unstable otherwise. Fluctuations in a are included through a stochastic 
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process ^(t), which we assume Gaussian, white and of small intensity k. The evolution of 
the system is now stochastic and is described by the joint probabihty density V{A, B; t) at 
time t. The reduction procedure starts by decomposing the joint density as 

P{A,B;t)=piB\A;t)PiA;t), (1) 

where p{B\A;t) is the conditional probability density. Close to threshold, the stochastic 
processes A and B are small (their intensity scales with some power of k) in such a way that 
characteristic values of-B/A~/t"^l, a>0. As will be shown in more detail below, this 
assumption also implies that the two processes evolve over different characteristic temporal 
scales, fact that is reminiscent of the separation of time scales present in the deterministic 
limit. As a consequence, the probability densities P{A;t) and p{B\A;t) can be separately 
obtained at different orders in n. The stationary density P{A) is then used to locate the 
effective threshold point in the stochastic case. Below threshold, P{A) is a delta function 
at y4 = 0, whereas above threshold there exists another normalizable solution that has some 
non vanishing moments. 

Van den Broeck et al. |^ introduced this reduction procedure to study the effect of 
additive noise on a pitchfork bifurcation. They derived an approximate expression for the 
stationary probability density near but below threshold. They showed that in the weak 
noise limit, the critical variable exhibits amplified non-Gaussian fluctuations and that the 
properties of the fast variable depend on the nonlinearity of the system under study. Their 
analysis, however, is difficult to extend to the region above threshold. We find that additive 
noise eliminates the separation in scales between the slow and fast variables, and that, as 
a consequence, the probability densities for A and B are in general quite broad. Hence the 
assumption that A/B 1 breaks down over significant portions of any particular trajectory, 
and the reduction procedure discussed is not reliable. 

In view of this limitation, the analysis presented here is restricted to equations involving 
multiplicative noise only. In this case, the separation in scales between the fast and slow 
variables is preserved well above onset. Knobloch and Wiesenfeld |^ had already addressed 
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the adiabatic elimination procedure in the multiphcative case by introducing one additional 
assumption: that fast variables are gaussianly distributed around the underlying determinis- 
tic center manifold. Our analysis extends theirs in that such an assumption is not necessary. 
In fact, we show that the fast variable does not always fluctuate around the manifold. 

We derive approximate expressions for the stationary probability densities p{B\A) and 
P{A) valid near threshold for the pitchfork and transcritical bifurcations. In both cases, the 
marginal density P{A) has to satisfy a normalizability condition that is used to determine 
the location of onset ac- In those cases in which etc 7^ 0, stochastic resonance between 
the fast variable B and the stochastic process ^{t) is responsible for the shift away from 
the deterministic threshold. This result generalizes earlier analyses of the normal form 
equation corresponding to a pitchfork bifurcation with a fluctuating control parameter [^,0 , 
in which coupling to fast variables was not considered. In agreement with our results below, 
the absence of such coupling leads to Oc = for any intensity of the fluctuating control 
parameter. 

The case of a pitchfork bifurcation with multiplicative noise is considered in Section 
For simplicity, the method is applied to the well-known Van der Pol-Duffing equation. In 
that example, the bifurcation point is shifted to ac > 0, while the fast variable B exhibits 
Gaussian fluctuations around B = 0. Our result for ac agrees with earlier work by Liicke 



10|, but disagrees with the work of Knobloch and Wiesenfeld M and of Seshadri et al. [11 



Section |T| considers the general case of a transcritical bifurcation. In this case, the fast 
variable can exhibit non Gaussian fluctuations and, in general, the mean of the distribution 
does not lie on the underlying deterministic center manifold. 

II. PITCHFORK BIFURCATION WITH MULTIPLICATIVE NOISE : THE VAN 

DER POL-DUFFING EQUATION 

In order to illustrate the reduction procedure in a model system that bifurcates super- 
critically, we consider the non-linear oscillator 
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also known as Van der Pol-Duffing oscillator The positive constants (3, a and b are of 
(9(1). At a = 0, Eq. (El) exhibits a supercritical pitchfork bifurcation between the two fixed 



point solutions x = (stable for a < 0) and x = ±ya/a (stable for a > 0). The last term 
in the right-hand side originates from a random component in the control parameter a. We 
limit our analysis to Gaussian, white noise satisfying (^(t)) = and {C,(t)C,(t')) = 2K6{t — t'), 
where ( ) denotes an ensemble average and k is the intensity of the noise. Motivated by 
the known center manifold reduction in the deterministic limit of k = 0, we perform the 
following linear change of variables A = x + x/ P and B = —x/P to yield 
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with c = a/ P,d = b — 3a/P,e = 26 — 3a/ P and / 

a/P a/P 



b — a/ p. The linear matrix L = 
has a zero eigenvalue at the deterministic bifurcation point of a = 0, 



-a/P -iP + a/P) 

with a second eigenvalue of 0{1). In the absence of noise, the variable B thus varies 
over a time scale which is much faster than the time scale of A. One then introduces the 
scalings a ~ C(e2), T = eH,A ~ 0(e) and B ~ 0{e^), with e « 1. Then, dB/dT ~ 
O(e^) ^ —{P + a/P)B + cA^, leading to the equation for the center manifold B„i{A) = 
{-aA/p + cA^)/{p + a/P) + 0{e^). Substituting this result into Eq. (|) gives the normal 
form equation for A. 

We now turn to the case k > 0, and keep the same change of variables under the 
assumption that the intensity of the noise is small: k ~ 0{e^). The exact Fokker-Planck 
equation associated with Eq. (H) is. 
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The first step in our analysis is to introduce scaled variables R = k,/ e"^ , A = A/ , B = B / 
and a = a/e^* in Eq. (Q). We choose a ~ in order to have a A ~ A'^ in the equation for 
A. In view of the deterministic result, we further assume that B/A ^ 1, and thus consider 
i < j. Eq. (H) now reads 
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Next, we introduce the decomposition V{A, B; t) = p{B\A] t)P{A] t) in Eq. (H) and integrate 
over B. Since [p{B\A;t)]^^^^ = [di^p{B\A;t)]i^^^^ = 0, all terms involving a derivative 
with respect to B integrate to zero, leaving 



dtPiA;t) = -e 
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Only the dominant contributions to the two terms on the right-hand side of Eq. (|^) were 
included. In order to obtain a dominant balance at (9(e^), we let i = 1. The marginal prob- 
ability density P{A] t) then evolves over a time scale T = t^t. By contrast, the conditional 
density p{B\A\ t) varies over times of 0{1), as seen from the equation 
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obtained by choosing j = 2 and restricting Eq. (|^) to 0{1). The separation of time 
scales central to the elimination procedure in deterministic systems is thus preserved in the 
stochastic case. The Langevin equation corresponding to Eq. (|^ is obtained by setting A 
to a constant in the original equation for B and dropping any term of 0{e^) or higher. 



The stationary solution to Eq. reads, in the original set of variables 
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It is a Gaussian distribution with zero mean and variance ct^(A) = Ky4^/(/5 + a/(3)(3'^. The 
fast variable -B thus fluctuates around S = and not around the center manifold Bm{A) (in 
contrast with the results of refs. and 0). In fact, Bm{A)/B ~ e'^/e^ ^ 1, thus indicating 
that terms proportional to aA and A'^ in the equation for B do not have any significant 
influence on its evolution. 

We note that the statistics of the fast variable are not generic but depend on the details 
of the system under consideration. For instance, if the equation for the fast variable is 
deterministic, the conditional density is a delta function on the center manifold [0. The 
procedure is then equivalent to replacing B in the equation for A by its value on the center 
manifold. Eq. (^) also fails if a term proportional to A"^ is present in the equation for B, in 
which case the Gaussian distribution is centered on the manifold -B^(A) ^ Const, x ~ 

The statistical properties of the critical variable A follow from Eq. (^. In particular, 
the stationary solution P{A) (or, equivalently, P{A)) to Eq. @ reads 



P{A) =^\A\^-^exp{-cf3^Ay2t^). (9) 



This density has nonzero moments and is normalizable (with J\f = [c/?^/2k] I /r(^ — i)) 



as 



long as ^ — 2 > — 1. This implies that, to C(e^), the bifurcation occurs at 

a, = K/p. (10) 

The bifurcation point is thus shifted to positive values of the reduced control parameter by 
an amount proportional to the noise intensity n. This result agrees with that of Liicke fi^ 
who used a perturbation analysis of the linear stability problem, but disagrees with earlier 



results due to Knobloch and Wiesenfeld ^ and Seshadri, West and Lindenberg |TT . 

We next compare our results (Eqs. (^ and (|^)) with a numerical integration of the 
original model equation (Eq. (0)). The numerical calculations were performed by using an 



explicit integration scheme, valid to first order in At |T3|, a step At = 0.005 and a bin size 
for the various probability densities AA = 0.01 and AB = 0.001. Initial conditions for x 
and X were chosen randomly from a uniform distribution in the interval [—0.5, 0.5]. Results 
from 100 independent runs were averaged, and within each run the various quantities were 
sampled every 1000 steps. To ensure the system had reached a stationary state, the first one 
million steps were discarded. Just above onset, the density Eq. exhibits a divergence 
at the origin (Fig. |l|A). At a = 2k,/ j3, this divergence transforms into a maximum (Fig. 
|1]B) which moves to higher values of A as the control parameter is further increased (Fig. 
IIP). All three figures, corresponding to the parameter values j3 = a = b = 1 and k = 0.01, 
show excellent agreement between the predictions of Eq. (|) and the stationary densities 
computed numerically. 

The average amplitude {\A\) was also computed for various values of a, and the results 
compared with the analytic result 
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r(a/3/2K- 1/2)' 

which follows from Eq. (|^). As shown in Fig. |I|D, agreement between the two data sets is 
once again excellent. As an additional test, we considered the statistics for the fast variable 
B. Combining Eqs. (|^) and (^ we find 



/oo i-oo / 2k 

^P(A)dAl Bp(BlA)dB = ^-^,^^^^m. (12) 

Analytic and simulation results are compared in Fig. ^ The quantity P{\B\) plotted in 
Figs. |]A, IpB and ^ is the density of probability of finding some value of independently 
of the value of A, i.e., P{\B\) = 2P{B) = Aj;^ P{A)p{B\A)dA. As before, both sets of 
results agree extremely well. 



III. TRANSCRITICAL BIFURCATION WITH MULTIPLICATIVE NOISE 

As a second illustration of the approach, we study the set of two equations 
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with a small and all the remaining coefficients of (9(1). In the deterministic limit, the 
variable B relaxes quickly to the center manifold Bfn{A) = dA'^/X, and the normal form 
equation is given by 

dA 
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which describes a transcritical bifurcation at a = 0. Following the procedure introduced 
above, we define the rescaled parameters R = n/e^ and a = a/e\ and rescaled variables 
A = A/e^ and B = B/e^, with i < j. The exact Fokker-Planck equation corresponding to 
Eq. (p!3D then reads 
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Integrating this equation over B gives 
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As in Section y, only the leading contributions to the right-hand side of Eq. (|T6|) were 
included. In order to have a dominant balance at 0(e^), we choose i = 2. Similarly, letting 
j = 3 in Eq. ( |15|) leads to the equation 



d,p{B\A-t) = -^[-XBpiB\A)] + ^[Rml.AMmi 
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valid to C(l). Eqs. (|^) and (|T^ admit the stationary solutions 
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with AT = (a/Rmj^) '^"^Ti ' ^Ti^ + ^^^^^^ + 1), and 



PiB\A) = J _\ exp 
y Z'KK.m2iA'^ 
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21^ 

respectively. The normalizability condition a/ztm^^ + 77112^21/^^1 — 1 > — 1 associated with 
Eq. ([TsD places the bifurcation point at 



etc = -^7771277721, (20) 

valid to 0{e'^). Predictions from Eq. (0) are compared with numerical estimates obtained 
through direct integration of Eq. ([T3|) in Fig. 0. The computations were performed with 
K, = 0.01 and all the parameters in the original equations except a set to one. The numerical 
and analytical estimates (represented by black dots and solid lines respectively) are virtually 
indistinguishable near onset. Significant differences do appear, however, as a increases. 
Results pertaining to the fast variable B are presented in Fig. ^. Again, analytic estimates 
of P{B) = Jf^"^ P{A)p{B\A)dA obtained by using Eqs. (|T|) and (|19|) compare well with 
their numerical counterparts near onset, but become increasingly inaccurate as a increases. 
In particular, the numerical results indicate that the density P{B) is slightly skewed and 
has a non-zero average. These properties are incompatible with a distribution such as Eq. 
(p!9|) which is even in B. 

We show next that it is in principle straightforward to systematically improve the accu- 
racy of the analytic calculation by going to higher orders in e. We seek a stationary solution 
of Eq. (|15D valid to one more order in e (C(e)). Setting dtV{A, B;t) = 0, i = 2, and j = 3 
in this equation and keeping terms up to C(e) gives, after some algebra, 

= -[-XB + edA"^ - e/t(777227772i + 3m2imu)A]P{A)p{B\A) 

+ [Kml^A^ + e(2/t7772 l 777 22^5)]P(^)-^^^^|^ + e(2/t777iir772l A^) 

oB 



(21) 
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In contrast with the calculation above, the equations for the conditional and marginal prob- 
abilities do not decouple. In order to solve Eq. (pTf ) for p{B\A), the derivatives dP{A)/dA 
and dp{B\A)/dA must be known to 0{1). 

We first determine dP{A)/dA by noting that the stationary solution to Eq. ([T6|) satisfies 
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We also assume that the conditional probability density p{B\A) is of the form 
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i.e., that the improved calculation simply adds corrections to the argument of the exponential 
in Eq. (|19]). Under that assumption, 
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Combining Eqs. (|21|), (|22|) and (0) gives, to C(e 
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The solution to that equation is an exponential, the argument of which can be expanded in 
the small quantity B/A. This yields the probability density 
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which is consistent with the assumption in Eq. 

Note that the presence of a cubic term in the exponential implies that the fast variable 
exhibits non Gaussian fluctuations. A divergence at either B ^ —oo or B ^ +oo also 
means that Eq. (^) is non-normalizable. In practice, however, one can compute an effective 
normalization constant by integrating Eq. (EBI) over some interval [B_,B^] at the limits of 
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which p{B±\A) <^ 1. Alternatively, higher order terms could be included in the Taylor series 
expansion leading to Eq. (p6D. For simplicity however, we let = rriii, in which case the 
coefficient of the cubic term vanishes and Eq. (^) simplifies to 

2am2i\ A'^ 



-e- [a + Kmi2m2i)A 

Amii 



B-e{d + 



mil ) A 
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Eq. (13) is a Gaussian distribution with mean {B) ^ = Bp{B\A)dB = e (d + ^ + 

^frnTT^"^ + Kmi2m2i)A different from the center manifold i?m(^), and variance {B^)a = 
RmliA^/X. 

As before, we determine the stationary properties of the slow variable A by setting 
dtV{A, B;t) = in Eq. ( [T5| ) and integrating over B. The resulting equation reads 

= ^ { [rttiI^A^ + e{2Rmiimi2A{B)A) + e^Rml^{B^)^ P{A)] 

a + Rml^ + Rmi2m2i)A - aA^ + e(2Kmiimi2 - - e^c{B'^)A P{A). (28) 

Inserting the expressions derived above for (-B)^ and {B^)a in that equation, and rearranging 
the various terms, we find 
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with g = —2m2i{a + Rmi2m2i) / Xrriii.r = {d + 2am2i / mn) / X and s = cRml-^^/X. An 
approximate solution can be found that, in the original variables, reads 

P{A) = AC'A'^ exp ( -^A^ - -^a] , (30) 
yKniii Kuiii J 
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and g' = — 2m2i(a + Kmi2m2i)/Amii. The dashed hnes in Figs. ^ and ^ are analytic 
estimates computed using Eqs. (|27| ) and (0). As expected, comparison with the numerical 
results shows a net improvement from our previous predictions (Eqs. (|T8|) and (p!9|)). 

We conclude our analysis by discussing briefly the mechanism by which the effective bifur- 
cation point differs from its deterministic location in the two cases studied above. Consider 
first the Van der Pol-Duffing oscillator, Eq. @, which we rewrite as 
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with a' = a/jS, mu = mi2 = and m2i = = — For simplicity, we only include in 
the deterministic part of Eq. (^) the terms which were found relevant in Section 0. Note 
that since A varies over a time scale T = e^t, the term proportional to in the equation 
governing its evolution can be averaged over the fast time scale. With the introduction of 
the scaled variables A = eA, B = e^B, a' = e^a', k = e^R and T = e^t, this gives 
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where we have approximated the temporal average of B^ by its ensemble average, and where 
(^(T)i^(T')) = 2e^R6{T — T'). By using the Furutsu-Novikov theorem M,|T5|, we find 



{B0 = {B){0 + e'R/^\=eRm2,A + 



(34) 



Therefore, the correlation of B^ itself evolves over the slow time scale. By combining Eqs. 
(p3|) and ( |3^ we obtain the effective normal form equation, 

dA 



dT 



dtA- A^ + e'^mnA^{T), 



(35) 



with a = a' + Km2imi2. This equation also leads to the Fokker-Planck equation (Eq. 
already derived in Section ||. The bifurcation point associated with Eq. (|35|) is located at 
a = Ul, i.e., at 



a = —Kmi2m2i, 
13 



(36) 



in agreement with our previous result (Eq. (|TOD)- Equation (^) is also identical to Eq. 



(pop derived in the case of a generic transcritical bifurcation. Hence, to first order in the 
intensity of the noise, the location of the bifurcation point is entirely determined from the 
stochastic part of the original set of equations and is therefore independent of the nature of 
the bifurcation. Corrections to Eq. (P^), however, depend on the details of the system under 
consideration. For instance, in the case of the transcritical bifurcation, the first correction to 
Eq. (|36D follows from the normalizability condition /i > — 1 associated with the probability 
density P{A) given in Eq. (^). 



ACKNOWLEDGMENTS 



This research was supported by the Microgravity Science and Applications Division of 
the NASA under contract No. NAGS- 1885, and also in part by the Supercomputer Com- 
putations Research Institute, which is partially funded by the U.S. Department of Energy, 
contract No. DE-FC05-85ER25000. 



14 



REFERENCES 

[1] J. Guckenheimer and P. Holmes, Nonlinear oscillations, dynamical systems, and bifur- 
cations of vector fields (Springer Verlag, New York, 1983). 

[2] R. Stratonovich, Topics in the Theory of Random Noise (Gordon and Breach, New 
York, 1967), Vol. II. 

[3] H. Haken, Synergetics: an introduction (Springer, Berlin, 1977). 

[4] C. van den Broeck, M. Malek Mansour, and F. Baras, J. Stat. Phys. 28, 557 (1982). 

[5] E. Knobloch and K. Wiesenfeld, J. Stat. Phys. 33, 611 (1983). 

[6] C. Elphick, M. Jeanneret, and E. Tirapegui, J. Stat. Phys. 48, 925 (1987). 

[7] F. Drolet and J. Viiials, Phys. Rev. E 56, 2649 (1997). 

[8] A. Schenzle and H. Brand, Phys. Rev. A 20, 1628 (1979). 

[9] R. Graham and A. Schenzle, Phys. Rev. A 25, 1731 (1982). 

[10] M. Liicke, in Noise in nonlinear dynamical systems, Vol. 2 of Theory of noise induced 
processes in special applications, edited by F. Moss and P. McClintock (Cambridge 
University Press, Cambridge, 1989). 

[11] V. Seshadri, B. West, and K. Lindenberg, Physica 107A, 219 (1981). 

[12] K. Wiesenfeld and E. Knobloch, Phys. Rev. A 26, 2946 (1982). 

[13] J. Sancho, M. S. Miguel, S. Katz, and J. Gunton, Phys. Rev. A 26, 1589 (1982). 

[14] K. Furutsu, J. Res. Nat. Bur. Standards 67D, 303 (1963). 

[15] E. Novikov, Sov. Phys. JETP 20, 1290 (1965). 



FIGURES 



0=0.015 



0=0.02 



30.0 



Numerical Simulations 
Predicted Density 




10.0 



• Numerical Simulations 
— Predicted Density 




I Numerical Simulations 
Predicted Average 




0.00 



0.05 



0.10 



|A| o 

FIG. 1. (A), (B) and (C): stationary probability densities as a function of the absolute value 
of A for the Van der Pol-Dufiing oscillator. Shown are the densities above onset for three different 
values of the control parameter a. (A), a = 0.015; (B), a = 0.02; and, (C), a = 0.04. We show in 
(D) the bifurcation diagram showing the average value (|^|) as a function of a. In all cases, the 
analytic results are represented by a solid line, whereas the symbols are the results of the numerical 
calculation. 
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FIG. 2. (A), (B) and (C): stationary probability densities as a function of the absolute value 
of B for the Van der Pol-DufHng oscillator. Shown are the densities above onset for three different 
values of the control parameter a. (A), a = 0.015; (B), a = 0.02; and, (C), a = 0.04. We show in 
(D) the bifurcation diagram showing the average value {\B\) as a function of a. In all cases, the 
analytic results are represented by a solid line, whereas the symbols are the results of the numerical 
calculation. 
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FIG. 3. (A), (B) and (C): stationary probability densities as a function of A for the transcritical 
bifurcation. Shown are the densities above onset for three different values of the control parameter 
a. (A), a = —0.005; (B), a = 0.01; and, (C), a = 0.04. We show in (D) the bifurcation diagram 
showing (A) as a function of a. In all cases, the analytic results to 0{1) and C(e^) are represented by 
solid and dashed lines respectively, whereas the symbols are the results of the numerical calculation. 
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FIG. 4. (A), (B) and (C): stationary probability densities as a function of B for the transcritical 
bifurcation. Shown are the densities above onset for three different values of the control parameter 
a. (A), a = —0.005; (B), a = 0.01; and, (C), a = 0.04. We show in (D) the bifurcation diagram 
showing {B) as a function of a. In all cases, the analytic results to 0(1) and 0(e) are represented by 
solid and dashed lines respectively, whereas the symbols are the results of the numerical calculation. 
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